from matplotlib.backends.backend_pdf import PdfPages
from astropy.visualization import MinMaxInterval, AsinhStretch, LinearStretch, \
        ImageNormalize, ManualInterval
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.ticker as ticker
import astropy.io.fits as fits
import astropy.stats
import numpy
import numpy as np
import glob
import os
import sys


#source = "IRAS15398"
source = 'V883 Ori'
v_source=4.25
print(source)


# Get a list of sources.



# A few definitions.

arcsec = 4.84813681e-06
c_l = 299792.458000

class Transform:
    def __init__(self, xmin, xmax, dx, fmt):
        self.xmin = xmin
        self.xmax = xmax
        self.dx = dx
        self.fmt = fmt

    def __call__(self, x, p):
        return self.fmt% ((x-(self.xmax-self.xmin)/2)*self.dx)

# Loop over field of view.
ticks_list = [numpy.array([-0.5,-0.3,0.01,0.3,0.5])]
velocity_range=[-4.5+v_source,4.5+v_source]
imagelist=['V883_Ori_SB_HDO-225.535GHz_robust_2.0.image.fits',
'V883_Ori_SB_HDO-241.558GHz_robust_2.0.image.fits',
'V883_Ori_SB_H218O-203.GHz-0.4kms_robust_2.0.image.fits']
keplerian_mask_list=[
'V883_Ori_SB_HDO-225.535GHz_robust_2.0.mask-hdo.fits',
'V883_Ori_SB_HDO-241.558GHz_robust_2.0.mask-hdo.fits',
'V883_Ori_SB_H218O-203.GHz-0.4kms_robust_2.0.mask-h218o.fits'
]
h218o_extra_mask_list=[
'V883_Ori_SB_H218O-203.GHz-0.4kms_robust_2.0.mask-ch3och3_blue.fits',
'V883_Ori_SB_H218O-203.GHz-0.4kms_robust_2.0.mask-ch3och3_red.fits',
'V883_Ori_SB_H218O-203.GHz-0.4kms_robust_2.0.mask-ch3och3_blue2.fits',
]
hdo241_extra_mask_list=[
'V883_Ori_SB_HDO-241.558GHz_robust_2.0.mask-ch3cho.fits'
#'V883_Ori_SB_HDO-241.558GHz_robust_2.0.mask-ch3cooh.fits'
]

hdo225_extra_mask_list=[
'V883_Ori_SB_HDO-225.535GHz_robust_2.0.mask-ch3ocho.fits'
]

rms_list=[1.4e-3,1.5e-3,1.55e-3]
levels_list=[np.arange(3,15,2)*rms_list[0],np.arange(3,15,2)*rms_list[1],np.arange(3,15,2)*rms_list[2]]
linenames=['HDO_225','HDO_241','H218O']
line_centers = [225.89672000,241.56155000,203.40752000]

continuum_image='V883_Ori_SB_continuum_235GHz_robust_0.5.image.tt0.fits'

# Loop through the sources and plot all the datasets for that source.

for image_name, line_center,line_name,mask_name in zip(imagelist, line_centers,linenames,keplerian_mask_list):
    # Now plot each of the images.
        # Load in the image
        print(line_name)
        image, header = fits.getdata("{0:s}".format(image_name), header=True)

        freq = numpy.arange(image.shape[0])*header["CDELT3"] + header["CRVAL3"]

        # Load the continuum image.
        data = fits.getdata(continuum_image)
        ndim_cont=len(data.shape)
     
        if ndim_cont == 4:
           cont = fits.getdata(continuum_image)[0,0]
        elif ndim_cont ==2:
           cont = fits.getdata(continuum_image) #[0,0]
        mask = fits.getdata(mask_name) #[0,0]
        # Get the center of the source(s).

        x0, y0 = 0., 0.

        # Check the number of rows and number of columns.
        velocity = (line_center - freq/1e9) / line_center * c_l

        velocities_to_plot=((velocity >velocity_range[0]) & (velocity < velocity_range[1])).nonzero()
        start_chan=velocities_to_plot[0][0]
        nrows, ncols = round(len(velocities_to_plot[0])**0.5)-1, round(len(velocities_to_plot[0])**0.5)+1
        #if line_name == 'HDO_225':
        #   nrows+=1
        #   #ncols+=1
        fontsize = 20
        rel_velocity=velocity -v_source
        v_sys_index=np.argmin(np.abs(rel_velocity))

        nearest_velocity_v_source=velocity[v_sys_index]
        
        print(rel_velocity)
        print(v_sys_index)
        print(nearest_velocity_v_source)

        for ticks in ticks_list:
            # Plot the image.

            N = image.shape[1]
            pixelsize = numpy.abs(header["CDELT2"])*numpy.pi/180. / arcsec

            xmin, xmax = int(round(N/2-x0/pixelsize+ticks[0]/pixelsize)), \
                    int(round(N/2-x0/pixelsize+ticks[-1]/pixelsize))
            ymin, ymax = int(round(N/2+y0/pixelsize+ticks[0]/pixelsize)), \
                    int(round(N/2+y0/pixelsize+ticks[-1]/pixelsize))

            npix = min(xmax - xmin, N)
            if xmin < 0:
                xmin, xmax = 0, npix
            if xmax > N:
                xmin, xmax = N - npix, N
            if ymin < 0:
                ymin, ymax = 0, npix
            if ymax > N:
                ymin, ymax = N - npix, N

            # Make a figure to put it in.

            fig, ax = plt.subplots(nrows=nrows, ncols=ncols, \
                    figsize=(2*ncols,2*nrows+0.5))

            # Get the normalization.

            vmin = -2*numpy.nanstd(image[0])
            vmax = max(10*numpy.nanstd(image[0]), numpy.nanmax(image))
            if line_name == 'H218O':
               vmax=vmax/2.0
            snr = vmax / astropy.stats.mad_std(image[0], ignore_nan=True)

            if snr > 100:
                norm = ImageNormalize(image, stretch=AsinhStretch(), \
                        interval=ManualInterval(vmin=vmin, vmax=vmax))
            else:
                norm = ImageNormalize(image, stretch=LinearStretch(), \
                        interval=ManualInterval(vmin=vmin, vmax=vmax))

            # Plot all of the channels.



            for i in range(nrows):
                for j in range(ncols):
                    ind = i*ncols + j +start_chan

                    # If we are beyond the limits of the frequency range, turn
                    # off the axis.

                    if ind >= image.shape[0]:
                        ax[i,j].set_axis_off()
                        continue

                    # Plot the data.

                    ax[i,j].imshow(image[ind,ymin:ymax,xmin:xmax], \
                            origin="lower", interpolation="bilinear", \
                            norm=norm, cmap="inferno")

                    # Contour the continuum data.

                    #ax[i,j].contour(cont[ymin:ymax,xmin:xmax], \
                    #        colors="white", levels=numpy.nanmax(cont)*\
                    #        (numpy.arange(5)+0.5)/5., linewidths=0.75, \
                    #        alpha=0.25)

                    #contour H218O image and mark channels with > 3sigma detection
                    if (line_name == 'H218O')        :
                       ax[i,j].contour(image[ind,ymin:ymax,xmin:xmax], \
                            colors="white", levels=levels_list[2], linewidths=1.1
                            )
                       if velocity[ind] > 5.3 and velocity[ind] < 6.7:
                          txt = ax[i,j].annotate(r"*".\
                            format(velocity[ind]), xy=(0.85,0.01), \
                            xycoords='axes fraction', fontsize=28, \
                            color="white",weight="bold")
                       if velocity[ind] > 1.9 and velocity[ind] < 3.5:
                          txt = ax[i,j].annotate(r"*".\
                            format(velocity[ind]), xy=(0.85,0.01), \
                            xycoords='axes fraction', fontsize=28, \
                            color="white",weight="bold")


                    #contour Keplerian Mask
                    if (line_name == 'HDO_241' and 8.2 > velocity[ind] > 0.1) or \
                       (line_name == 'HDO_225' and 8.2 > velocity[ind] > 0.1) or \
                       (line_name == 'H218O' and 8.2 > velocity[ind] > 0.1)        :
                       ax[i,j].contour(mask[ind,ymin:ymax,xmin:xmax], \
                            colors="white", levels=numpy.array([1.0]), linewidths=1.5
                            )

                    #contour Keplerian Mask
                    if (line_name == 'H218O')        :
                       mask_red = fits.getdata(h218o_extra_mask_list[1]) #[0,0]
                       mask_blue = fits.getdata(h218o_extra_mask_list[0]) #[0,0]
                       mask_blue2 = fits.getdata(h218o_extra_mask_list[2]) #[0,0]
                       ax[i,j].contour(mask_red[ind,ymin:ymax,xmin:xmax], \
                            colors="red", levels=numpy.array([1.0]), linewidths=1.5
                            )
                       ax[i,j].contour(mask_blue[ind,ymin:ymax,xmin:xmax], \
                            colors="blue", levels=numpy.array([1.0]), linewidths=1.5
                            )
                       ax[i,j].contour(mask_blue2[ind,ymin:ymax,xmin:xmax], \
                            colors="blue", levels=numpy.array([1.0]), linewidths=1.5
                            )
                    if (line_name == 'HDO_241')        :
                       mask_blue = fits.getdata(hdo241_extra_mask_list[0]) #[0,0]
                       ax[i,j].contour(mask_blue[ind,ymin:ymax,xmin:xmax], \
                            colors="blue", levels=numpy.array([1.0]), linewidths=1.5
                            )
                    if (line_name == 'HDO_225')        :
                       mask_blue = fits.getdata(hdo225_extra_mask_list[0]) #[0,0]
                       ax[i,j].contour(mask_blue[ind,ymin:ymax,xmin:xmax], \
                            colors="blue", levels=numpy.array([1.0]), linewidths=1.5
                            )
                   # Add cross to mark the continuum position

                    txt = ax[i,j].annotate(r"+", xy=(0.455,0.465), \
                            xycoords='axes fraction', fontsize=fontsize, \
                            color="white")

                    # Add labels to the x-axis.

                    transform = ticker.FuncFormatter(Transform(xmin, xmax, \
                            pixelsize, '%.1f'))
                    transform2 = ticker.FuncFormatter(Transform(xmax, xmin, \
                            pixelsize, '%.1f'))

                    ax[i,j].set_xticks((ticks[1:-1]-ticks[0])/pixelsize)
                    ax[i,j].set_yticks((ticks[1:-1]-ticks[0])/pixelsize)
                    ax[i,j].get_xaxis().set_major_formatter(transform)
                    ax[i,j].get_yaxis().set_major_formatter(transform)
                    #a work around since transform function doesn't want to plot reverse
                    ax[i,j].set_xticklabels(['0.3','0.0','-0.3'])
                    if i == nrows-1:
                        ax[i,j].set_xlabel('$\Delta$R.A. (arcsec)', \
                                fontsize=fontsize/1.75, labelpad=8)
                    else:
                        ax[i,j].set_xticklabels([])
                    if j == 0:
                        ax[i,j].set_ylabel('$\Delta$Dec. (arcsec)', \
                                fontsize=fontsize/1.75, labelpad=12)
                    else:
                        ax[i,j].set_yticklabels([])

                    ax[i,j].tick_params(axis='both', direction='in', \
                            labelsize=fontsize/1.75, color="white")


                    # Make the axes white as well.

                    for side in ["left","right","top","bottom"]:
                        ax[i,j].spines[side].set_color("white")

                    # Mark V_sys channel
                    #if velocity[ind] == nearest_velocity_v_source:
                    if ind == 10+start_chan or ind == 11+start_chan:
                       txt = ax[i,j].annotate(r"$\bigstar$",\
                             xy=(0.8,0.8), \
                            xycoords='axes fraction', fontsize=fontsize*1.5, \
                            color="white")
                    # Add the velocity to the image.
                    txt = ax[i,j].annotate(r"$v={0:3.2f}$".\
                            format(velocity[ind]), xy=(0.01,0.85), \
                            xycoords='axes fraction', fontsize=fontsize, \
                            color="white")
                      
                    # Show the beam.

                    bmaj = header["BMAJ"]/abs(header["CDELT1"])
                    bmin = header["BMIN"]/abs(header["CDELT1"])
                    bpa = header["BPA"]

                    xy = ((xmax - xmin)*0.1, (ymax - ymin)*0.1)

                    ax[i,j].add_artist(patches.Ellipse(xy=xy, width=bmaj, \
                            height=bmin, angle=(bpa+90), facecolor="white", \
                            edgecolor="black"))

            # Add a title to the figure.
            if line_name=='HDO_225':
               titlename='HDO (3$_{1,2}$-2$_{2,1}$) (225.896 GHz)'
            if line_name=='HDO_241':
               titlename='HDO (2$_{1,1}$-2$_{1,2}$) (241.561 GHz)'
            if line_name=='H218O':
               titlename='H$_2$$^{18}$O (3$_{1,3}$-2$_{2,0}$) (203.407 GHz)'

            if nrows*ncols > 30:
               fig.suptitle("{0:s}".format(titlename), fontsize=2*fontsize)
            else:
               fig.suptitle("{0:s}".format(titlename), fontsize=fontsize)

            # Adjust the spacing.

            plt.tight_layout(pad=0, rect=(0.01,0.01,0.99,0.95))

            # Save the figure.

            plt.savefig(line_name+'.pdf',dpi=200.0)

            # Clear the figure.

            plt.close(fig)

